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PERFECT PLASTICITY WITH DAMAGE AND HEALING AT 
SMALL STRAINS, ITS MODELLING, ANALYSIS, AND 
COMPUTER IMPLEMENTATION 

TOMAS ROUBi'CEK* * * § tt AND JAN VALDMANi§ 

Abstract. The quasistatic, Prandtl-Reuss perfect plasticity at small strains is combined with 
a gradient, reversible (i.e. admitting healing) damage which influences both the elastic moduli 
and the yield stress. Existence of weak solutions of the resulted system of variational inequalities 
is proved by a suitable fractional-step discretisation in time with guaranteed numericalstability 
and convergence. After finite-element approximation, this scheme is computationally implemented 
and illustrative 2-dimensional simulations are performed. The model allows e.g. for application in 
geophysical modelling of re-occurring rupture of lithospheric faults. Resulted incremental problems 
are solved in MATLAB by quasi-Newton method to resolve elastoplasticity component of the 
solution while damage component is obtained by solution of a quadratic programming problem. 
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1. Introduction. There is a vast amount of literature about plasticity and 
about damage separately, both in mathematics and in civil or mechanical engineer¬ 
ing. Much less literature addresses various combination of plasticity and damage, 
cf. e.g. [SUSHUjlinilllJIlSIST]. In engineering, this is usually called ductile damage, 
cf. e.g. 1 18lf28Ll30tl35j . Also a lot of geophysical models combine reversible damage 
(called rather ageing) with some sort of plasticity (often modelled as not entirely 
independent of damage, however), cf. e.g. [521 . 

The goal of this article is to devise a model that would allow for 

• modelling of thin plastic shear bands surrounded by wider damage zones (as 
typically occurs in geophysical modeling of lithospheric faults with very narrow 
core) with possible healing of damage (as considered in geophysical modeling 
to allow re-occurring damaging), and simultaneously 

• rigorous proof of existence of weak solutions of the resulted system of varia¬ 
tional inequalities proved by a suitable fractional-step discretisation in time 
with guaranteed numerical stability and convergence, and 

• efficient numerical implementation of the time-discrete model. 

We depart from the standard linearized, associative, rate-independent plasticity at 
small strain as presented e.g. in [25]. Simultaneously, we use also a rather standard 
scalar (i.e. isotropic) damage as introduced by L.M. Kachanov in late 60ieth and 
presented e.g. in m, considered here however as rate dependent and reversible in 
the sense that a possible healing is allowed. To avoid serious mathematical and 
computation difficulties, we have in mind primarily an incomplete damage through 
a higher-order damage-independent term, although the standard elastic tensor can 
allow for a complete damage, cf. H and C = C(C) below. An important aspect 
of the model is that not only the conservative part but also the dissipative part is 
subjected to damage, i.e. not only the elastic moduli but also the yield stress will be 
considered as damageable. This relatively simple and lucid mechanism will however 
lead to a possibly very complex response of the model. 


* Mathematical Institute, Charles University, Sokolovska 83, CZ-186 75 Praha 8, Czech Rep. 

^Institute of Thermomechanics, Czech Acad. Sci., Dolejskova 5, 182 00 Praha 8, Czech Rep. 

$ Institute of Information Theory and Automation, Czech Academy of Sciences, Pod 
vodarenskou vezf 4, CZ-18208 Praha 8, Czech Republic. 

§ Institute of Mathematics and Biomathematics, Faculty of Science, University of South Bo¬ 
hemia, Branisovska 31, CZ-370 05 Ceske Budejovice, Czech Republic. 


1 



2 


T.Roubfcek & J.Valdman 


To make the model accessible to analysis, we work within the setting of small 
strains, and we also take into account surface-energy effects by including in the free 
energy a term dependent on the gradient of the total strain. This is also known as a 
concept of so-called second-grade nonsimple materials, cf. e.g. [40[l50l . alternatively 
also referred as the concept of hyper- or couple-stresses (42[[5lj ; for reasons we use 
it here cf. Remark 12.51 below. 

In view of applications we have in mind, we suppress any hardening effects and 
thus we consider the Prandtl-Reuss elastic /perfectly plastic model; in fact, consider¬ 
ing kinematic or isotropic hardening would make a lot of aspects even much easier. 
A plastic yield stress dependent on damage is in some variants used in the Cam-Clay 
model, cf. e.g. [T21IHI1HE], or in the Perzyna model with damage, cf. [5T[, and also 
in (2l[3lf9li~RTj . Let us also point out that damage with healing without plasticity (as 
sometimes considered in mathematical literature) would have only very limitted ap¬ 
plication because damaged material typically can undergo substantial deformation 
and the healing should not be performed towards the original configuration. 

We confine on the isothermal variant of the model. In contrast to [48] , we 
consider rate-independent plasticity without any gradient, so that concentration of 
plastic and total strains and development of sharp shear bands is possible. Also, 
related to this concentration, both plastihcation and damage are driven by the 
elastic stress (which is still well controlled) rather than the total strain (which may 
concentrate); for plasticity itself, see also [47] . 

The presented model has potential application in geophysical modelling of re¬ 
occurring rupture of lithospheric faults or of nucleation of new faults. A narrow 
so-called core of the fault can be modelled by the perfect plasticity while and a 
relatively wide damage zone around it can arise by the gradient-damage model. 
After a combination with inertial effects (and possibly a visco-elastic rheology e.g. 
of Jeffreys type), this model involves seismic waves and can serve for earthquake 
simulations where these waves are emitted during fast rupture, cf. Remarks 12.31 and 
12.41 below for some modifications of the presented model towards these applications. 
Another possible modification, going beyond the scope of this paper however, might 
use the structure of the stored energy similar to what is used in a phenomenological 
models for polycrystalline shape-memory alloys where our damage variable is in a 
position of temperature and plastic strain is a transformation strain subjected to 
some additional constraints, see e.g. [HI Example 5.15]. 

The plan of the paper is as follows: In Section [2] we formulate the model and 
cast a suitable definition of the weak solution, and pronounce a basic existence result 
which is proved later in Sections [3] by a constructive time discretisation method. 
A further finite-element discretisation is then outlined. This allows for computer 
implementation of the model presented in Section [4[ whose efficiency and some 
physical aspects eventually demonstrated on in Section [5] an illustrative example 
with geophysical motivation. 

2. The model, its weak formulation, and existence result. Hereafter, 
we suppose that the damageable elasto-plastic body occupies a bounded smooth 
domain H C R d , d = 2 or 3. We denote by n the outward unit normal to dCl. We 
further suppose that the boundary of H splits as 


dn ■.= r = r D u r N , 


with r D and T N open subsets in the relative topology of d f2, disjoint one from each 
other, each of them with a smooth ((d^l)-dimensional) boundary, and covering <9f 1 
up to (<J—l)-dimensional zero measure. Considering T > 0 a fixed time horizon, we 
set 

Q:=(0,T)xn, E:=(0,T)xT, S D := (0, T)xT D , E N := (0, T)xT N . 
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Further, and R dx v d will denote the set of symmetric or symmetric trace-free 

(= deviatoric) (dxd)-matrices, respectively. For readers’ convenience, let us sum¬ 
marize the basic notation used in what follows: 


d = 2, 3 dimension of the problem, 

Ri X v d == {-4 G R; tr A = 0}, 
u : Q — > R d displacement, 

7r : Q — y Rde/ plastic strain, 

( : Q — > [0, 1] damage variable, 
a : R —¥ R + damage-dissipation potential, 
b : [0,1] — > R stored energy of damage, 
e e i elastic strain, e e i = e(u)—Tr, 
e = e(u) = |Vu T + 

total small-strain tensor, 

C : [0,1] —)• R 3 elasticity tensor 
dependent on £, 


f| hyperstress (3rd-order) tensor 
H a (small) hyperelasticity tensor, 

5 = a Y (-)Si [0,1] =t Rl x v d , 

with B\ the unit ball in Rj X v d , 
a Y : [0,1] —> R + plastic yield stress 
dependent on 

g : Q —► R d applied bulk force, 

Wn '■ E D —> R d prescribed time-dependent 
boundary displacement, 

/ : E n —¥ R d applied traction force, 
k > 0 scale coefficient 

of the gradient of damage. 


Table 1. Summary of the basic notation used thorough the paper. 

The state is formed by the triple q := (u,ir, £)• Considering still a (small but fixed) 
regularizing parameter e > 0, the governing equation/inclusions read as: 


(2.1a) div(C(C)e e i — div (i) + g = 0 (momentum equilibrium) 

with f ) = HVe e i and e e \ = e(u)— tt, 

(2.1b) dSg^(n) 9 dev(C(C)e e i — div t)), (plastic flow rule) 

(2.1c) da{() + ic'(C)e e i : e e i 

- Kdiv((l+ £ |VCr~ 2 )VC) +iV [0 ,i](C) 9 b'(Q, (damage flow rule) 

with 5s the indicator function to S and 5g its convex conjugate. Here, [C(C)e]y 
and [UVe]^^ mean and . n _-| 1HI r .j mn e%n ■ respectively. 

We employed two regularizing terms with a regularizing tensor El and a regular¬ 
izing parameter e > 0 with an exponent to be assumed suitably big, namely r > d. 
This regularization facilitates analytical well-posedness of the problem and, because 
the gradient-damage term degenerates at VC = 0, its influence is presumably small 
if e is small and VC not too large. Moreover, H in (1 2. tail prevents a complete dam¬ 
age at least when we assume C(C) positive semidefinite. Actually, (|2.1b|) represents 
rather the thermodynamical-force balance governing damage evolution while the 
corresponding flow rule is written rather in the (equivalent) form 

tt e N S (q (dev(C(C)e e i - div ft)) 

with N the set-valued normal-cone mapping to the convex set indicated. An anal¬ 
ogous remark applies to (I2.1cl) . 

A remarkable attribute of this model is a damage-dependent yield-stress do¬ 
main S = 5(C). Typically, developing damage makes S smaller and vice versa, i.e. 
S(-) : [0,1] Rdev^ nondecreasing with respect to the ordering of subsets by 
inclusion. Likewise, typically also &(•) and C(-) are nondecreasing, the later one 
with respect to the Lowner’s ordering, i.e. C(zi) — 0 ( 22 ) is positive semi-definite 
for zi > z 2 . Rate-dependency of damage evolution prevents nonphysically too-early 
damaging/plastification and, due to the driving force b'( C), also allows simply for re¬ 
verse damage evolution (a so-called healing) by using a convex function a : R — > R + 
in (I2.1cl) having naturally its minimum at 0. The microstructural interpretation of 
b is a stored energy related with microcracks/microvoids arising by damage, reflect¬ 
ing the fact that any surface in the bulk bears some extra energy. Minimization of 
this energy naturally leads to a tendency for healing of these material defects. Of 
course, EH> is to be completed by appropriate boundary conditions for (I2.1h .c). 
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e.g. 



(2.2a) 

u = W D 

on r D 

(2.2b) 

(C(C)e e i - divf))-n- diVg(fjn) = / 

on r N 

(2.2c) 

VC-n = 0 and f):(n®n) = 0 

on r 


with n denoting the unit outward normal to £1. Moreover, div s is the surface- 
divergence operator, which may be introduced as follows [22]: given a vector field 
v : T —► R d , we extend it to a neighborhood of T, and we let its surface gradient 
(valued in M. dxd ) be defined as V s v = P s Vw, where P s = I — ft <S> n is the projector 
on the tangent space of T; we then let the surface divergence of v be the scalar 
field div s i; = P s : V s m = tr(P s VvP s ). Given a tensor field A : T —> R dxd , we let 
div s A : r —»• be the unique vector Held such that div s (A T a) = a-div s A for all 

constant vector fields a : T —> Furthermore, the symbols “ • ” and “ : ” denote a 

contraction between the one or two indices, respectively. Later, we will use also “ j” 
for a contraction between three indices. Thus, componentwise, the second condition 

in (12.2b|) reads as J2j,k= 1 ^ijknjUk = 0. _ 

Of course, an inhomogeneous variant of (12.21)1) or some mixed Dirichlct/Neumann 
conditions in the normal/tangent conditions could be considered with straightfor¬ 
ward modifications of the following text. We will consider an initial-value problem 
for (12.11) (12.21) by asking for 

(2.3) u(0) = tt( 0) = 7T 0 , and C(0) = Co¬ 

in fact, as u does not occur in m, u o is rather formal and will essentially be 
determined by 7r 0 and Co via (12.14hD below. 

The system m with the boundary conditions (12.21) has, in its weak formula¬ 
tion, the structure of an abstract Biot equation (or here rather inclusion): 

(2.4) d^(q;q) + d£{t,q) B 0 

with suitable time-dependent stored-energy functional § and the state-dependent 
(pseudo)potential of dissipative forces Equally, one can write (12.41) as a gener¬ 
alized gradient flow 

(2.5) q e dfM* (q; -d£{t, q)) 

where C | -t &*{q', C) denotes the conjugate functional to n-> ^(g; v). 

The perfect-plasticity model itself received considerable attention already a long 
time ago, see e.g. in [5l lTTlH4ll26ll351l44] . The peculiarity is that the displacement 
no longer lives in the conventional Sobolev H 1 -space but rather in the space of 
functions with bounded deformations introduced by Suquet [53], defined as 

(2.6) BD (n-,R d ) := {ueL\n;R d ); e{u) £ Meas^; K^)}, 

where Meas(fi) = denotes the space of Borel measures on the closure of fb 

The other notation we will use is rather standard: beside the standard notation for 
the Lebesgue L p -space we already used in (12.61) for p = 1, we further use W k,p for 
Sobolev space whose fc-th derivatives are in L p -spaces, the abbreviation H k = W k ’ 2 , 
and L p (0,T; X) for Bochner spaces of Bochner-measurable mappings (0,T) —> X 
with X a Banach space. Also, W k,p ( 0, T; X) denotes the Banach space of mappings 
from L p { 0, T;X) whose fc-th distributional derivative in time is also in L p ( 0, T ; X). 
Further, C([0,T];X) and C' wea k([0, T]; X) will denote the Banach space of continu¬ 
ous and weakly continuous mappings [0, T) —> X , respectively. Moreover, we denote 
by BV([0, T];X) the Banach space of the mappings [0, T] —> X that have a bounded 
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variation on [0, T], and by B([0, T]; X) the space of Bochner measurable, everywhere 
defined, and bounded mappings [0,T] —> X. 

After considering smooth time-dependent Dirichlet boundary conditions w D on 
£ d which allows for an extension onto Q , let us denote it by u D , such that 

(2.7a) (C(C)e(u D ) — div f) D )-n — div s ((j D iu) = 0 on r N , 

(2.7b) f) D :(n ® ft) = 0 with l) D = HVe(w D ) on T 

for any admissible C, and making a substitution of u + u D instead of u into EH 
(12.21) . we arrive to the problem with time-constant (even homogeneous) Dirichlet 
boundary conditions. More specifically, 

(2.8a) e e \ in (12.1bl) replaces by e e i = e{u+u D )— tt, and 

(2.8b) w D in (12.2al) replaces by 0. 

The state space is then the Banach space 

(2.9a) U:= { (u, tt, C) £ BD(D; R d ) xMeas(D; M(j e x v d )x W 1 ’ r (D); 

e(u)— 7 r € i 7 1 ( D ; Kg y x I ((), u 0 ndS+ tt = 0 on T D }, 


where a 0 b means the symmetrized tensorial product -j(a ® b + b ® a), and the 
functionals governing the problem (12.41) leading to (12.11) (12.21) with the substitution 


(12.81) are: 

(2.9b) 


(2.9c) 


£{t,u,w, C) := < 


[ ^C(C)(e(u+u D (i))-7r) : (e(u+u D (t))- 7r) 

In * 

+-HV(e(u+« D (f))-7r): V(e(u+u D (t))—7 t) 

-b{0 ~ 9 {t)-u + k (||VC| 2 +^|VCr) da; 

— f f(t)-udS if [0,1] a.e. on D, 

ar N 

oo otherwise, 


@(£’,*,0 : = / [^( c ) M ] ( da: ) + / a ( C ) d £, 

Jn Jn 


where S$/q denotes the conjugate to the indicator function 8s(q to the convex set 
S(C) and where the first integral in (12.9c|) is an integral of a Borel measure; counting 
the assumption (12.1 dfl) below, this measure is cr Y (C)|7r| with |7r| the total variation 
of tt. The norm on U is 


|| (u,TT, C)|| u IMIl^R 1 *) + ll e ( M )llMeas(fi;R^) 

+ 11 11 Meas(n ; Rj^ d ) + \\ e ( u )~HI H i (fijR**?) + IICII W 1 > r (n) • 

We can now state the weak formulation of the initial-boundary-value problem 
EH EH- As for the plastic part, we use the concept of the so-called energetic 
solution devised by Mielke and Theil i59], cf. also El [S3, based on the energy 
(in)equality and the so-called stability and further employed in the viscous con¬ 
text in |05] with the stability condition modified to a semi-stability, cf. (12.Hal) 
below. Another feature of the following definition is that we rely on a regular¬ 
ity of the damage C so that div((l+e|VC| T '~ 2 )^ 7 C) is in duality with ( and thus, 
in fact, the damage flow rule (12.lcD holds even a.e. Q. Actually, we do not need 
such regularity for the definition itself because the usual weak formulation of (|2.1cl) . 

which would involve (not well-controlled) VC resulted from usage of Green’s for¬ 
mula, could be still treated by applying a by-part integration in time to get rid off 
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the term ((l+e|VC| r_2 )VC) • V£. Rather, this regularity is essential for the energy 
conservation. 

Definition 2.1 (Weak solution). The triple (u,ir,£) with 
(2.10a) u £ B([0,T];BD(fi;R d )), 

(2.10b) tt e B([0, T]; Meas(; K^* v d )) fl BV([0, T]; Meas(Q; M dx v d )), 

(2.10c) C € B([0, T]; W lir (fi)) n H 1 { 0, T; L 2 (D)) n C([0, T] xO) 

such that also 


(2.10d) e e i = e(u+n D ) — 7T £ B([0,T];if 1 (D;R dxd )) and 

( 2 . 10 e) div((l+e|VCr“ 2 )VC) € L 2 (Q) 

is called a weak solution to the initial-boundary-value problem (EH) EH with the 
substitution (12.81) if: 

(i) the semi-stability 


(2.11a) £(t, u{t),n(t), C(t)) < &{t, u, tt, C(i)) + tt - n(t), 0) 

holds for all t £ [0,T] and for all £ BD(D; R d ) xMeas(D; Kj* v d ) with 

uQndS+Tr = 0 on T D and with e(u)—ir € H 1 (f2;R d ^), 

(ii) the variational inequality 


(2.11b) 



i<C'(C)e el : eei 


-Kdiv((l+e|VCr~ 2 )VC) 

_ fc'(C) + €) {v - c) dxdt > 


a(f) da: df, 


holds for all v£ L 2 (Q) and some f£L 2 (Q) such that ^£ a.e. onQ, 

(Hi) the energy equality 


(2.11c) £(T,u(T),n(T)X(T))+ [ [5* s(<) (n)](dxdt) + [ a(Q dxdt 

l[0,T]xfI JQ 

= #(0,wo,7ro,Co) + f d t £(t,u(t),Tr(t)X(t))dt. 

Jo 

holds with a : R —> R being the single-valued, continuous function defined by 
a(z ) := zda(z). 

(iv) and also the initial conditions EH hold. 


Let us note that, counting cancellation of some terms in d?(t,u(t),n(t),£(t)) — 
$(t, u, tt, C(f)), the semi-stability (12.llal) means that 

( 2 - 12 ) ^C(C(t))(e(u(t)+2u D (t))-7r(t)(t)) : (e(u(t))-n(t)) 

+ -HV(e(u(t)+2u D (t))—7r(t)(t)) : V(e(u(f))—7r(t)) da; 

< J ^C(C(t))(e(u+2u D (t))-7r) : (e(u)-Tr) 

+ (e(u+2u D (t))—n) :V(e(u)-n) dx + J_ [5s K(t)) (7r-7r(t))] (da:). 

The last integral (12.121) is not a Lebesgue integral but an integral according the 
measure S* Bi (fk— 7r(t)). Due to the special ansatz (12.14 fl) below, this integral will the 
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total variation \tt— 7r(f)|, namely f s ~ } cr Y (C(f))|7T—7r(t)|(dx). Similarly, the integral 
on the left-hand side of (12.llell equals Jj Q T ] x q (J Y (^)|7r|(da;dt). Further note that, 
although traces of functions from BD(f2;R d ) are in L 1 (r;R d ), one has to be aware 
of jumps that can occur at the boundary, i.e. the measure e(u) may concentrate on 
the boundary T. Thus, the classical boundary condition u = 0 on r D arising by the 
additive shift (I2.8bl) is replaced by the more involved relation u 0 ndS + n = 0 on 
r D in (12.flail . This relation has to be understood as an equality of measures on T D : 

V measurable A C T D : uQ ndS = dn = it {A). 

J A J A 

The relation simply means that any jump of u on the boundary has to be due to 
a localized plastic deformation. Cf. m for analytical details. Eventually, let us 
comment the last term in (12. 1 lcl) which, in view of (I2.9bl) . involves the expression 

(2.13) d t £(t,u, 7r, C) = f C(£)(e(w+w D (f))-7r) : e(u D (t)) 

Jn 

+ iV(e(M+w D (t))-7r): Ve(n D (f)) — g(t)-udx — / f(t)-udS. 

J r N 

Let us collect the assumptions on the data and on the loading we will rely on, 
some of them being already mentioned above: 

(2.14a) f 1 C R d bounded C 2 -domain, T D has a (d— 2) dimensional C' 2 -boundary, 
(2.14b) a : R —> R convex, smooth on R\{0}, a(0) = 0, and 
3e > 0 Vz€R : e|^:| 2 < a{z) < (l+|. 2 | 2 )/e, 

(2.14c) b : [0,1] —> R continuously differentiable, non-decreasing, concave, 

(2.14d) C : [0,1] —> R dxdxdxd continuously differentiable, positive-semidefinite, 

V b ji l = 1; * * • , d I &ijld (•) = Cjikl (*) = &klij (*), 

VeeRg^ : C(-)e:e : [0,1] —¥ R non-decreasing, convex, 

3C d (C), c s (C) : C(C)e : e = C D (C)deve : deve + c s (C)(tr e) 2 , 

(2.14e) H positive definite, HI ijki = H jiM = H k uj, 

3H D , H s : HVe : Ve = H D Vdeve : Vdeve + H S V tr e • Vtr e, 

(2.14f) S(C) = a Y : [0,1] —>■ (0,oo) continuous nondecreasing, 

with Bi C Rje V d a unit ball, 

(2.14g) w D eW 1A (0,T-H 3/2 (T D ;R d )) and 3 u D € W 1A {0, T; i7 2 (fl; R d )) 

satisfying (12.711 and n D |r D = *di 

9 £ W 1 - 1 (0,T;L 1 (n ; R d )), / € W 1 ' 1 (0,T-,L 1 {T t ,;R d )), 

3 cr SL : [0, T] —> L 2 (VL\ Rg^m) 3 a > 0 : a SL n = g on[0,T]xT N and 
divcr SL + / = 0 and |devcr SL | < <r Y (0) — a on[0,T]xfl, 

(2.14h) (n 0 ,t r 0 , Co) £ BD(fi;R d )xMeas(fi;R^ x v d )xW 1 ’ r (0), 

0 < Co < 1 a.e. on Q., and 
V(u,7f) £ BD(n ; R d )xMeas(fi;R^ e x v d ), 

e(u)— tt Si7 1 (f2; Rf^), u 0 nd5+ tt = 0 on T D : 

^(0, Mo, 7T 0 , Co) < ^(0, u, TT, Co) + ^(Co; 0, 7T—7T 0 ), 

(2.14i) k > 0, £ > 0, r > d. 


The smoothness assumption (12.14a.ll and the “elastic” invariance of the orthogonal 
subspaces of deviatoric and volumetric components (I2.14l l.e) copy the assumptions 
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used in D33 for perfect plasticity in simple materials without damage in a variant 
with spatially varying yield stress as in mmm- The stress cr SL in the condition 
P-14g[ ) qualifies the loading be / and g in such a way so that the infinite sliding 
of some parts of body is excluded; this is a usual requirement called a safe-load 
condition, connected to perfect plasticity, here adopted to the situation that the 
yield stress a Y may vary with damage similarly as in [15 [ Remark 2.9]. It should be 
also remarked that this safe-load condition works similarly for nonsimple materials. 
Further note that (12.1 Till) represents in particular the semi-stability of the initial 
condition and makes, with other assumption, the energy conservation (12.11cl) pos¬ 
sible. Note also that (12.14bll ensures that a used (12.11cl) is single-valued although a 
itself may be set-valued at 0. In (12.14H) . one can easily consider a bit more general 
situation when B\ would be convex, closed, and 0 £ int B \. 

The main analytical result justifying rigorously the model (12.11) (12.31) is: 

Theorem 2.2. Under the assumptions (|2.14|) . at least one weak solution to the 
initial-boundary-value problem m (usd according to Definition \2.1\ does exist. 

We will prove this existence result in Section [3] by a constructive time discreti¬ 
sation method, cf. Lemma 13.11 with Proposition 13.31 which later in Sections 0] and 
© allows for efficient computer implementation of the model. The uniqueness of 
the solution however hardly can be expected. 

Remark 2.3 (The dynamical model). During fast rupture, inertial effects may 
be not negligible and even sometimes an important aspect of the model. Then, 
(12.1 all augments by the inertial force gu with g > 0 denoting the mass density as 

(2.15) gu — diver = g with er = C(<C)e e i — div t). 

Relying on that the inertial term gu is controlled in the space L 2 (0, T; H 2 (Cl\ R 3 )*) 
D CweakQO, T\; L 2 (Ul] R 3 )) or actually even in a slightly better space counting that 
dever £ L°°(Q;]Rj^), the weak formulation of (12.151) arising by double by-part 
integration in time should accompany (12.111) with $ augmented by the inertial en¬ 
ergy f Q ||u| 2 da; but with (12.1 lal) holding only a.e. on [0,T] and (12.1 lcl) only as an 
inequality. The functional in (13.5al) then augments by gr~ 2 \u — 2 u k ~ l + u k ~ 2 1 2 /2. 
Actually, it seems a matter of a physically-explainable fact that some difficulties the 
energy conservation occurs probably due to integration of elastic waves with non- 
linearly responding shear bands even if a Kelvin-Voigt-type visco-elastic rheology 
would be involved, cf. also m Remark 6]. In this dynamical case, the fast dam¬ 
age phases and subsequent fast plastic slips, called (tectonic) earthquakes, typically 
emit elastic (seismic) waves. However, although some justification on theoretical 
level, the computational modelling requires fine special techniques to suppress e.g. 
parasitic numerical attenuation and the direct combination of elastic waves with 
the inelastic processes is difficult. 

Remark 2.4 (A non-Hookean model). The concept of nonsimple materials 
allows an important generalization that S'ft,-,^,-) is not quadratic and even non- 
convex. More specifically, instead of the coercive term (e e i, C) l “t <C(£)e e i:e e i = 
+ n(C)h as used also here in (14.11) below, [33] proposed 

(2.16) (e e i,C) ^A(C)/i +^(0-^2 with h = tre e i, I 2 = |e e i| 2 - 

The elastic stress is then (A(C)— 7 (C)v^ 2 )t re ei + (2/r(£)— 7 (C)-fi/V^ 2 ) e eb while the 
driving stress for damage is crdam = \X'(QI 2 + //(CK 2 — 'y'iOh/Vh and can now 
be positive even without the contribution of the 6-term. Such a model is widely 
used in geophysics where it is believed to be responsible for instability of heavily 
damaged rocks and leads to healing even without the 6-term used in our model, 
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but where it is used without the H-term and thus without any rigorous justification 
of such models, cf. e.g. [251151] and references there. To preserve coercivity of the 
model due to boundary conditions and the H-term, one can think about a certain 
softening under very large strain by replacing 2-homogeneous form (12.1611 by an 
energy with only a linear growth 


(2.17) 


(e e i, C) 


A(C)J 1 2 +2 M (C)/ 2 -2 7 (Q/iv / l2 


with e > 0 presumably small. A certain conceptual inconsistency remains in 
damage-dependence of C but not of H, although H is assumed to be only small 
in applications. Note that (I3.5al) then represents a coercive but non-convex min¬ 
imization problem and one should seek a global minimizer to ensure (!3.9aP . The 
nonsimple-material concept allows for a simple modification of the convergence proof 
in semistability and in the damage flow by compactness: more specifically, the bino¬ 
mial trick in (13.171) is applied only to the dissipation and the H-terms, while (13.181) 
is even simpler because C'(£)e e i is now bounded in L°°(Q;W lxd ). 

Remark 2.5 (A simple-material model). Considering H = 0 would bring vari¬ 
ous difficulties. In particular, the L 2 (Q)-estimate of the driving force iC'(C)e e i:e e i, 
which would need a regularity of e e i that however does not seem available for plas¬ 
ticity models without hardening, would become problematic. Note that the higher 
integrability of e e i®e e i will be used e.g. in (13.181) and in (13.211) too. One should note 
that the alternative idea to consider a nonlinear damage independent contribution 
to the stress of the type +£|e e i| 2 e e i would not allow to use the binomial trick in 
the Step 3 in the proof of Proposition 13.31 below, while the strong convergence of 
e e i seems also not obvious to prove. A certain possibility might be in considering a 
visco-elastic Kelvin-Voigt model with the stress H>(C)e e i-|-C(£, e e i) with a nonlinear, 
monotone C(£, •) having at most the growth |C(£, e e i)| < C{ 1 + leeil 1 / 2 ) so that 
f 0 d(C((, te e i) d t can still be estimated in L 2 (Q) due to the D-term which can even 
depend on f as in 158 ] , 

3. The discretisation, its stability and convergence. To implement the 
initial-boundary-value problem (12.IP (12.3P computationally, we need to make a time 
and space discretisation. 

Let us first make only a time discretisation with, for notational simplicity, a 
constant time step r > 0. As the inertial effects are not considered and thus the 
system is only lst-order in time, the dependence of r > 0 on the time levels is easy 
to consider for numerical analysis and to implement (as actually used in Section [I] 
below). 

As $ is convex in terms of (it, 7r) and separately in £ too, and also as £% additively 
splits (ft, 7r) from £, the natural fractional-step strategy leading to an efficient and 
numerically stable semi-implicit formula follows this splitting (it, 7r) from £. More 


specifically, it reads 

as 


(3.1a) 

div(c(C T fc - 


O 

II 

+ 

1 


with 

4, ,T = 

e{u k +u D (kT))-TT k , f^=HVe^ r , g k T := g{kr) 

(3.1b) 

N S{ C?- 1 ) ( 

'K k — 'K k 

T 

-) 9dev(c(C^)e e fc 1>T -div(^), 

(3.1c) 

/C k -C k 
da(^ T W 

-) + 

1c'(#) 4 iT = 4,r 




K div((l+ £ |VC T fc r- 2 )VC T fe ) +N m (<: k ) 9 b\C), 
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together with the corresponding boundary conditions 

(3.2a) u k = 0 on T d , 

(3.2b) (C(d _1 )eg lT - divf)£) • n- div s (t) k n) = f k on T N with f k := /(fcr), 
(3.2c) V( k -n = 0 and t)£:(n®n) = 0 on T, 

to be solved first for (u k ,n k ) from (13.lh .bY (13.2h .bl and then for ( k from (13.1cl) 
()3.2cl) recursively for k = 1 ,...,T/r. Both these boundary-value problems have 
potentials and thus leads to minimization problems. Moreover, as C and —b' are 
nondecreasing (again with respect to the Lowner’s ordering) and a is convex as 
assumed in (12.1411 , both these boundary-value problems leads to convex variational 
problems, cf. (13.51) below. 

Let us define the piecewise affine interpolant u T by 
(3.3a) u T (t):=- ——— ^-u k + — —- u k_1 for t £ [(fc— l)r, kr] 

T T 


with k = 0, Besides, we define also the left-continuous piecewise constant 

interpolant u T and the right-continuous piecewise constant interpolant u T by 


(3.3b) u T (t) := u k for t £ ((fc— l)r, kr] , k = l, ...,T/r, 

(3.3c) u T (t) := u^ 1 for t £ [(k— l)r, kr) , k = l,...,T/r. 

Similarly, we define also n T , W T , 7r r , £ T , ( r , ffr, etc. 

Lemma 3.1 (Existence and stability of discrete solutions). The recursive boun¬ 
dary-value problem (ED m has a weak solution (u k , n k , £^) with u k £ BD(f2; R d ), 
7 x k £ W lir (fi), and Cr e Meas(n;R^ x v d ) with e k Xr = e{u k )~n k £ H 1 { f2;Rf y x ^) sat¬ 
isfying the a-priori estimates 


(3.4a) 

(3.4b) 

(3.4c) 

(3.4d) 

(3.4e) 


^ T IL“(0,T;BD(n;R d )) — ^ i 

7r ' r llz,“>(0,T;Meas(n;R^ e >< v <i )) nBVQO.TJjLitnjR^)) — ^ 

e (u T )— 7 rT|| ioO ( 0 )T;irl ( n . R dxd)) < c, 

< ’ T IL'»(o,T;W 1 >'-(n))nff 1 (o,T;L 2 (n)) — ^ 
div((i+e|vc r r- 2 )vc r )|| i 2 (Q) <c. 


Proof. The existence of weak solutions to EED can be justified by the di¬ 
rect method when realizing the variational structure: the boundary-value problem 
(l3.lh .bV (13.2h .bl represents a minimization problem 

{ Minimize (u, 7 r) i-a S{kr 1 u, 7r, C^ _1 ) + &{( k ~ 1 ]ix—Tx k ~ 1 ,0) 

subject to u £ BD(fi; R d ), 7r £ Meas(fi; R“[ x v d ), 

e(u)— 7r£ H 1 (0; Rg y x ,)(), u 0 ndS 1 -!- 7r = 0 on T D , 

while the boundary-value problem (I3.1cl) (I3.2cl) represents a minimization problem 


(3.5b) 


) Minimize 

g(kT,u k ,Tr k X) + T&(( k x ;0, 

1 

H ?r 

1 

1 subject to 

C £ W lir (n), 0 < C < 1 on Q, 



whose solutions do exist by coercivity, convexity, and lower semicontinuity argu¬ 
ments. Here the safe-load qualification (|2.14gj) of / and g is to be used. 
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Further, we test (13.11) respectively by u*—u* _1 , n k — 7r* _1 , and £* —£* _1 . Re¬ 
lying on the convexity of £(kr,-,-,( k ~ 1 ) and of £{kr, u k , 7r*, •), we obtain the esti¬ 
mates 


it*" 1 ,^- 1 ,^- 1 ), 


(3.6a) £{kT,u k ,n k ,Cr *) + [a- Y (C 1 )|7r*—7r* 1 |(dx) < £ (hr, 

Ja 

(3.6b) <?(fcr, «*,**,£) + [ atf-g-^dx < ^(fcr.u*,^,#" 1 ) 

J n 

with a from (I2.11el) . By summing these estimates, we can enjoy the cancellation of 
the terms <f(fcr, it*, 7r*, £* _1 ) hr (13.6aD and (I3.6b|) . and we thus obtain 


(3.7) ^(fcr,^, 7 r*,c?)+^(cr 1 ;<-<" 1 .c;-cr 1 ) < ^u*- 1 ,**- 1 ,# -1 ) 

nkr 

= £((k-l)T,u k T -\n k T -\C~ 1 )+ / W,«‘- 1 ) 7rJ- 1 ,C^- 1 )dt 

J (k— l)r 


with the dissipation rate ^ defined as 

(3.8) ^(C;tt,C):= [ cr Y (C)|7r|(da;) + [ a{() dx with a(C) = C#a(C). 

Jn Jn 

By summing (13.71) over k we enjoy a “telescopic” cancellation effect. Realizing (12.131) 
and p,14g[ ), by the discrete Gronwall inequality, we obtain (13.4b .-dl. 

Having estimated da(( T ) + !C'((()eei,T : e e i, T — £/(C r ) as a bounded set in L 2 (Q) 
uniformly with respect to r > 0, we can estimate also div((l+e|VC*|’" _2 )V((*) in 
the same space. For this, we test (13.1cl) by —div((l+e|VC*r _2 )VC*). Here, the 
important ingredient is, written rather formally, the following estimate 

/ iv [01] (c*)(-div((i+ e |vc*r _2 )vd)) dx 

J n f 

= - / a<5 [0jl] (c*)(div((i +£ |vdr- 2 )vc T fc ))dx 
Jn 

= [ v(^ [0 .i](d))-(i+eivc T fc r- 2 )vddx 
J n 

= / 5 2 <5 [0il] (C T fc ) • Vd-(l+e|VC?r 2 )Vd dx > 0 

Jn 

which is due to the positive-semidefiniteness of the (generalized) Jacobian d 2 d^ 0 -n 
of the convex function <5^ 0 ^ and which is to be proved rigorously by a mollification 
of & [ 0 , 1 ], cf. 03 Lemma l] for technical details. Thus we obtain (I3.4cl) . □ 

Lemma 3.2 (Discrete analog of (12.111) 1. With the notation (13.31) and e e i, T = 
e(u T +u DyT ) — 7T t , the discrete solution obtained by the recursive scheme (ED ED 
satisfies: 

(3.9a) £(t,u T (t),W T {t),C r {t)) < £(t,u,n,C T (t)) +&(C T (t)-,Tv - W T (t),0) 

for all t £ [0, T\ and all as in (12.11al) . and 


( 3 . 9 b) f a(v) + (ic'(CJe e i, T : e el , T - k div((l+ e |VC r r 2 )VC T ) 

Q , _ _ v • f • 

— b\( T ) + f T ) (v - ( T ) dxdt > / a(C T )dxdt 
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holds for all v £ L 2 (Q) and for some £ T £ L 2 (Q) such that £ ^V[o,i] (C T ) a - e ■ on 
Q, and eventually the energy (im)balance holds: 


(3.9c) £(T,Ur(T),n T (T),Cr(T))+ [ &(t:n T ,( T )dt 

Jo 

d t £(t,u T (t),n T (t),C T {t))dt 

with the overall dissipation rate 2% from (13.81) . Moreover, the a-priori estimate 
holds: 


< <^(0,uo,7r o ,Co) + [ 
Jo 


(3.10) 


kr\ 


< a 

L 2 (Q) - L ■ 


Proof. The boundary-value problem (|3.1h .b)- (|3.2h .b) represents a minimiza¬ 
tion problem (|3.5aD which can be tested by (u* -1 ,and, by using a triangle 
inequality facilitated by the 1-homogeneity of •, C), we obtain (13.9a|) ; actually, 
this is a standard argument in the theory of rate-independent processes [36j|371[39]. 

In the case of the boundary-value problem (13.1cl) (13.2cl) , the variational inequal¬ 
ity (I3.9bl) represents just the conventional weak formulation of the minimization 
problem (I3.5bl) summed for all time levels. Then, (13.9cl) follows by summing (13.71) 
for k = 1,..., T/t. 

Eventually, the estimate (13.101) follows by comparison from the inclusion £ 

fo'(Cr) ~ l < ^' , (C T )®ei,r : 5 e i )T + k div((l+e| V ( T | 7 ' _2 ) VC T ) — da(( T ) and by the already 
obtained estimates. □ 

Proposition 3.3 (Convergence). Let the assumptions (12.141) be satisfied and 
the approximate solution (u T ,W T ,£ T ,£ T ) be obtained by the recursive scheme «o 
m- Then there is a subsequence and (u,tt, £, £,) such that 

(3.11a) u T (t) —»• u(t) weakly* in BD(fi;R d ), 

(3.11b) W T (t) —> n(t) weakly* in Meas(fi; Rj* v d ), 

(3.11c) e e i, r (i) = e(u T (t))—W T (t) —> e(u(f))—7r(f) = e e i(t) weakly in U 1 (fl; R^J), 
(3.lid) ( T (t) —> f(t) and ( r (t) —> £(f) weakly in W 1,r {M) 

holding for any t£ [0,T], and further also 

(3.lie) f f strongly in L°°(Q), and 

(3.Ilf) £ r —> f weakly in L 2 {Q) 

with from Lemma \3.2l Moreover, any obtained by such a way is a weak 

solution according Definition \2.1\ with £ in (I2.11bl) taken from (I3.11fl) . 

Proof. For clarity of exposition, we divide the proof into five particular steps. 

Step 1: Selection of a converging subsequence. By Banach’s selection principle, 
we select a weakly* converging subsequence with respect to the norms from the 
estimates (ECU) and (13.101) : namely, for some u, i r, C, and f we have 

(3.12a) u T ->• u weakly* in L°°(0, T; BD(fi; R d )), 

(3.12b) 7f T ^ 7T weakly* in L°°(0, T; Meas(0; R(j e x v d )) n BV([0, T]; L 1 ^; R^)), 
(3.12c) e e i, r = e(u T )— tt t —> e e i = e(u )—k weakly* in L°°(0, T; H 1 ^- R^J)), 
(3.12d) Cr -► C weakly* in L°°(0, T; lT 1 ’ r (fl)) D H\0, T; L 2 (n)), 

(3.12e) div((l+e|VC T r" 2 )VC T ) ->■ div((l+ e |VCr“ 2 )VC) weakly in L 2 (Q), 
(3.12f) £ T —► £ weakly in L 2 (Q)\ 
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actually, (I3.12cl) uses also the maximal monotonicity of the involved nonlinear op¬ 
erator. Moreover, by the BV-estimates and the Helly’s selection principle, we can 
also count with (I3.11bl) and £ T (t) —>■ ((t) weakly in L 2 (ft), and then by the a-priori 
W 1,r -estimate (I3.4dh also both the first and the second convergence in f|3. 1 Idfl : both 
limits in (I3.11dll are actually the same because the limit ( is continuous in time into 
L 2 (Vl) due to the the a-priori iJ 1 -estimate <|3.4df> . 

By the compact embedding kF 1,r (fi) <s C(f2) and by the Arzela-Ascoli modifi¬ 
cation of the Aubin-Lions theorem, cf. [46], Lemma 7.10], we have the compact em¬ 
bedding CweadM; W 1 - 1 ^)) n JT 1 (0,T;L 2 (n)) <§ C([0,T]-,C{ti)) = C{Q). Thus, 
from the estimate (I3.4dl) . we obtain Cr —t C i n C(Q). Further, we have 


(3.13) 


CrLo=( 0 ,T;i=(n)) = SU P / \C T (^ X ) - Cr{t,x)\dx 

< [ ( sup |C (t,x) - Cr(i,x)| 2 )da; 

Jn v o <t<T J 



>fc -i 


2 dx < 



rk-i\ 


'dx 



2 

d x = T 



dxdt. 


Then, using the Gagliardo-Nirenberg inequality ||z|| ioo(n) < C , e ||z||| j2(n) ||x;||^ 1 % (n) 
for some small 0 < e < 1 depending on r > d, we can interpolate (13.1311 . i.e. 

IIC T — Cr||i=°(o,T;L 2 (n)) < V^ llCr|[ l 2 (Q) , with ||C T — Cr||(o,T ; ivc r (n)) < C to obtain 
||C T - CrII l°°(Q) -t o. Thus (13.11c|) is proved. 


Step 2: Energy inequality. The convergence (13.121) allows already for passage in the 
limit in the inequality (13.9cl) by lower semicontinuity in the left-hand side and by 
continuity in the right-hand side of (13.9cl) . 

The limit passage in £{T, u t {T), i t t (T), ( t (T)) is by the convexity of £(T, 
and the compactness in (, while for d t £{t,u T {t), 7r T (t), £^(t)) dt we use the con¬ 
tinuity of d t £(t,■,■,•) from (12.131) and the Lebesgue theorem; more in detail, we use 
the assumptions Q2.14gj ) and the weak convergence (13.11cl) . 

The only remaining (and nontrivial) term is the dissipation t^-term. Let us 
note that, as the discrete flow rule N s ^ ^(n T ) 9 dev(C(C T )e e i,r — div \) Y ) as well as 

the dissipation rate <x Y {C, T )\'h T \ uses C T and not just ( T , we needed to prove (13.11cl) 
in Step 1. Therefore, we have at disposal the estimate 


(3.14) 


K CT y(C t ) O'y (C) ) I 111 Meas (Q) - ^<t y ||C t C|| L oo(Q) ||^|| M eas(Q) 


with £ the modulus of Lipschitz continuity of cr Y on [0,1], cf. the assumption 
(I2.14fl) . Then, using also ( T —> f in C(Q) already proved, we obtain 


(3.15) liminf / ; n T , ( T ) dt = liminf / <t y (£ ) I 7 t t I (dxdt) 

T ^° Jo ~ T T ^° Jq ~ t 

= lirn J (c y (C t ) — <t y (£)) |7 Tt-| (dxdt) + liminf J cr Y (£)|7T r | (dzdt) 


> 0 + / CT Y (C)|d|(da;dt); 

Jo 


for the used weak* lower semicontinuity of n 1-9 JqCt y (C)|tt|( drcdi) we refer e.g. 

to mm- 
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Step 3: Limit passage in the semi-stability (|3.9ap towards (|2.11af) . For any (u,tt) 
used in (I3.9al) . we have to find at least one so-called mutual recovery sequence 
{(rt T , 7r r )} r >o in the sense that 

limsup £(t, U T , 7 Tr, ( (t)) + (t); 7f T -7f T (t), 0) - £{t, Ur{t),W T (t), C (*)) 

r-S-0 

< £{t, u, tt, CM) + r (t); 7 T- 7 T (t), 0 ) - £(t, u{t),TT(t), c (t)). 

We choose 

(3.16) u T = u T (t) + u — u(t) and 7r r = 7f r (t) + 7r — 7r(i). 

Then, by using the cancellation and the binomial formula of the type a 2 — b 2 = 
(a+b)(a—b) here in the form like Ce:e — Ce:e = C(e+e):(e—e) and HVe-Ve— 

HVe -Ve = HV(e+e) ;V(e—e), cf. (12.121) . and by making the substitution (13.161) . we 
have 


(3.17) lim £(t, u T , 7r r , C{t)) + &(C(t); 7r T -7r T {t), 0) - £{t, u T (t), tt t ( t), C(t)) 

T —>-0 — T — T — T 

= lim (^J ic(C T (t))(e(u T +u T (t)+2u D (t))-7f T -7f T (t)) 

: (e(Ur-Ur{t))-TTr+Tfr(t)) - g(t)-(u T -U T (t)) 

+ -HV {e(u T +u T (t)+2u D (t))— tt t — n T (t)) 

: V {e(u T — U T (t))~ TTr+TTr (t)) dx 

+ f_[(T v (C {t))\n T -7f T {t)\(dx) - [ f(t)-(u T -Ur{t))ds] 

J n J r N J 

= lim ^ J ic(C T (t)){e( : u T +u T (t)+2'u D (t))-7f r -7f r (t)) 

: (e(u—u(t))—Tr+Tr(t)) — g(t)-(u—u(t)) 

+ ^M'V[e(u T +Ur(t)+2u D (t))—Trr~Tfr(t)) : V {e(u—u(t))— 7r-|-7r(t)) dx 

+ [ VY(( T (.t))\n-ir(t)\{dx)) - [ f(t)-(u-u(t))dS 
Jn J J r N 

= / -C(C(t))(e(u+5(t)+2u D (t))-7r-7f(t)) : (e(u—u(t))—n+ir(t)) 

Jn 2. 

+ —MV(e(u+u(t)+2u D (t))— it— 7f(f)) : V(e(u—u(t))—n+7r(t)) dx 

+ A ^(Ol^r—7r(^)| (d^) - f g(t)-(u-u(t))dx - [ f(t)-(u-u(t)) dS 

J £1 J J In 

= £{t, U, vf, CM) + ^(CMi 7f—7r(t) } 0) - (t, u(t), 


Note that we used also Q' Y (M M)I 71 ~—^Ml ~ cr Y (C) |?r—7r(i)| in Meas(fl) due to the 
continuity assumption 12. 1411) on er Y and due to the convergence CM —► CM in 
C(f2) which follows from the second estimates in (I3.11dl) and the compact embed¬ 
ding w 1 - r (n) c c(n). 

Limit passage in the damage flow rule (13.9b!) towards (12. 1 lb!) . We need 
to prove that e e i, T —> e e i strongly in L 2 (Q; Rf^). To this goal, we first realize 
that Ve e i, r (t) —> Ve e i(t) weakly in L 2 (fl;R dxdxd ) as pronounced in (13.11clk here 
we use the uniqueness of stresses (counting the already selected subsequence (13.121) 
and its limit), cf. the arguments in [TTJ Thm.5.9] or also in [351 Sect.4.2.3] for 
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simple materials without damage. Here, using also absolute continuity valid due to 
viscosity in damage flow rule we obtain 


(3.18) V( e “-e™)> + <C(<)(ei I) -e< 1 2) ),e< 1 I) -e' 1 2 >>) 




< max |C'(z)| IICll 7 - 2 / 0 ) I 
o<z<i 'un"L 2 (n)i 


„(i) 


a ( 2 ) II 2 

'cl IIl 4 (Q;I 


Note that, for H = 0 and C' = 0, it reduces to the simple inequality (a^ — 
®ei ^ — ^e\) — 0 use d ' n Here, we should integrate (13.181) over [0,i], 

use positive-definiteness of H and C(-), and eventually Gronwall’s inequality, which 
works here certainly even for d < 4 for which the embedding U 2 (B) C 
holds. By this way, we obtain = e^. Thus, using the compact embed¬ 
ding, we also know that e e i, T (t) —> e e i (t) strongly in L 6 _e (H;if d < 3. 
Then, by the uniform bounds in time and by Lebesgue’s theorem used e.g. to 
t ||e e i, T (t) - e e i(t)|| L1 (n-s dxd )i we can see that e e i jT —> e e \ strongly even in 
L e ~ e {n;R^)) with each small e > 0. 

Then the only difficult remaining terms are k Jq div((l+e| V£ T | r-2 )VC T )Ci- dxdt 

and Jq dxdt because so far we know only the weak convergence of of 

div((l+e|VC r | r_2 )VC r ), and of £ T in L 2 (Q). We indeed cannot expect the limit, 
but we can proceed the following estimate: 


(3.19) limsup f div((l+e|VC T r 2 )VC t )Ct dxdt 

r-s-0 Jq 

= — liminf / (l+£r|VC T r _2 )VC r -VC T dxdt 

T ^° Jq 

< limsup / i|VCo| 2 + -|VCor-^|VC.(T)| 2 --|VC T (T)| r dx 

r—>o Jq 2 r 2 1 r 

< f J|VCo| 2 + f |VC o r-J|VC(T)| 2 --|VC(T)| r dx 

Jq 2 r 2 1 r 


= / div((l+e|VCr~ 2 )VC)Cda;dt 


JQ 

where we used (13. 1 lcll) at t = T and where the last equality relies on the regularity 
property div((l+e|VC| r_2 )VC) £ L 2 {Q) and can be proved either by a mollification 
in space 1411 Formula (3.69)1 and or in time by a time-difference technique 1211 
Formula (2.15)]. 

The convergence in the inclusion £ T £ ./V^p] (C r ) is easy due to the maximal 
monotonicity of lV[ 0 .i](-) and the convergences (13. 1 lfl) and ( T —>■ ( strongly in L 2 (Q) 
which can be proved by a generalized version of the Aubin-Lions theorem, cf. j461 
Corollary 7.9], or here even in L°°(Q) was proved as in Step 1. Having proved 
£ £ A r [ 01 ]((C), we can also see that 


(3.20) limsup [ £ T (-< T ) dxdt = limsup ( / <5 [0 i](Co) dz - / <Wi](CrC r )) dz 
t^o Jq r—>o \Jq Jq 


< / d [0 . 1] (Co)dz- / d [01] (C(T))dz = / £(-<() dzdt, 
Jq Jq Jq 


which is needed for the limit passage in (I3.9bl) : in fact, even the limit and the 
equality hold in (13.201) . 












16 


T.Roubfcek & J.Valdman 


Step 5: Energy equality. We test (12.1cl) which holds a.e. on Q by (. This test is 
legal as all terms in (12.1 cl) as well as £ are in L 2 (Q). We again use the last equality 
in (13.191) . Moreover, as £ € <9<5[ 0 q(£), we have fqfC dxdt = f Q 5[ 0) i] (£(T)) — 
^[o,i](C(0))dx = 0 — 0 = 0. We thus obtain 

(3.21) [ ^|vC(T)| 2 + -|VC(T)| r -6(C(T))dx 

do 2 1 r 

+ f ^C'(C)e e i : e e i+d(C)dxdt = / ^|VCo| 2 + —|VC 0 |' - KCo) dx. 

Jq 2 do 21 r 

Furthermore, we test formally (12.lal) by u and (|2.1bl) by fi. The rigorous cal¬ 
culations uses the approximation of the Stieltjes-type integral by Riemann sums 
and semistability, cf. m Formulas (76)—(82)] which adapts technique developed 
in the theory of rate-independent processes mini. Here, as C is not constant, 

we will still see the term (^C^^eepeeOC which results by the formal substitution 

C(C)e e i:e e i = Jj§C(C)e e i:e e i - (dC'(C)e e i:e e i)C; note that C(C)e e i:e e i is not well de¬ 
fined since e e i is not well controlled. Thus we obtain 

(3.22) [ ic(C(T))e e i(T):e el (T) + i]HVe el (T) i Ve e i(T) dx 

do 1 1 

+ [ _CT Y (C)|7r|(dxdt) = [ (ic'(C)e e i:e e i)cdxdt 
d[o,T]xO Jq 7 

+ / ic(Co)e e i(0):e e i(0) + ilHIVeei(O): Ve e i(0)dx. 
do 1 1 

Summing (13.211) and (13.221) then gives the energy balance (12.1 lcl) . □ 

Further, to implement the model computationally, we need to make a spatial 
discretisation of the time-discrete scheme m (Ei- To this goal, we use the 
lowest-order conformal finite-element method (FEM). In view of the used regular¬ 
ity (I3.4cp . the straightforward discretisation therefore employs P2-elements for u 
and C and Pl-elements for n. Rigorously speaking, due to the assumed smooth¬ 
ness (I2.14a| . one should consider FEM on a nonpolyhedral, curved domain. The 
minimization problems (13.51) are then to be restricted on the corresponding finite¬ 
dimensional subspaces, and the solution thus obtained is denoted by u* h , and 
0 h , with h > 0 denoting the mesh size. By this way, we obtain also the piecewise 
constant and affine interpolants in time, denoted by u T h and u T h , 7f T h and n T h, 
and eventually ( Th and ( T h■ Also, £ rh can be obtained analogously as before in 
Lemma 13.21 

Proposition 3.4 (Convergence of the FEM discretisation). Let (|2. 141) be sat¬ 
isfied, and the P2-FEM for u and C, and Pl-FEM for n is applied with h > 0 the 
mesh size. Then: 

(i) the a-priori estimates (13.41) and (13.101) hold when modified for u T h, 7r T h, ( T h, 
and £ Th with C independent of t > 0 and now of h > 0, too. 

(ii) Moreover, these fully discrete solutions converge (in terms of subsequences) 
in the mode as (13.111) towards weak solutions according Definition \2.1\ when 
simultaneously r —> 0 and h 0. 

The modification of the proof of this joint convergence of time-and-space dis¬ 
cretisation is rather routine, the explicit construction of the mutual recovery se¬ 
quence (13.161) taking additionally a finite-element approximation like in [5] , namely 
u T h = u Th (t) + n^ 2) (u - u{t)) and it Th = fr Th {t) + TL^ (?r - tt (f)) with nj^ and n^ 2) 
denoting a projector onto the PI- and P2 FE-spaces, respectively; we omit details 
about this modification. 
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Remark 3.5 ( Damage discretised by P 1-elements). The damage flow rule 
(12.lei) itself suggests to use only P 1-elements for £ which are, naturally, more easy 
to implement than the P2-elements used in Proposition 13.41 Then however (I3.4cl) 
cannot be expected for the FEM approximation of ( and also a direct P-1 FEM 
analog of (I3.9bl) cannot hold. Instead of (I3.9bl) . we have 

( 3 - 23 ) (a{v) + (^ C '(C Tft )ee i,r : e e i,r - &'( Crh) + trh) (v - Crh) 

+ K ((l+elV( rh r 2 )V( Th )-V(v-( Th )) dxdt> J a(( Th )dxdt 

for any v valued in the finite-dimensional Pl-FE subspace. Yet, the sequence 
{VCrh} r>o,h>o cannot be expected bounded. Thus, for the limit passage, instead of 
(13.231) one should rather use the discrete by-part integration (summation) in time 
like we did in (13.191) . leading to 

( 3 - 24 ) J C (aty) + (C Th )eelT ■ eel,. ~ 6'( C Th ) + C T h) ('V - Crh) 

+ « ((l+ e |VC Th | r - 2 )VC Th ) • dxdt + ^ ||VCo| 2 + y |VCor dx 
> [ a(Crh)dxdt+ [ ?-\V( Th (T)\ 2 + — \\7( Th (T)\ r dx 

Jq Jn 2 r 

which holds for any v valued in the P 1-finite-element space. Now, however, we do 
not have the estimates (13.4cjl and (13.101) . Anyhow, the limit passage seems possible 
by using the strategy proposed by Colli and Visintin [8]> cf. also [46] Sect. 11.1.2], 
allowing for the stored energy S taking values +oo but relying on boundedness of 
as indeed our situation. The convergence is, of course, in a weaker mode than 
(13.111) . Only after this limit passage, we can prove the regularity (12.Kiel) and go 
back to the weak formulation (I2.11bl) by using also the arguments which we use for 
the last equality in (13.191) . 

4. Implementation of the fully discrete model. The implementation of 
the model addressed in Proposition [T4] is rather cumbersome because of high-order 
FEM involved. Therefore we dare make few shortcuts: Pl-elements are used for 
damage ( according to Remark 13.51 Moreover, the (anyhow usual small and even 
not reliably known) hyperelasticity moduli are neglected, i.e. H = 0 and then small- 
strain tensor gradients Ve(u) are not involved. Consequently, only Pl-elements can 
be used for displacement u and PO-elements for plastic strain. Only the case d = 2 
is treated, so the previous analytical part have required r > 2 and we dare make 
another (indeed small) shortcut by considering r = 2 (and therefore by putting 
e = 0 the damage-gradient term in (I2.9bl) become quadratic). 

The material is assumed isotropic with properties linearly dependent on dam¬ 
age. The isotropic elasticity tensor is assumed as 

(4.1) ^dijkl(C) ■— [(^1 ^o)C Y Ao \dijSkl “f~ [(/T 1 po)C Y po \(dik$jl Y SilSjk) 

where Ai, fi i and Ao, po are two sets of Lame parameters satisfying 

Ai > A 0 > 0, pi > o Y 0. 

Here, <5 denotes the Kronecker symbol. This choice implies that the clastic-moduli 
tensor satisfies (I2.14dl) and it is even positive-definite-valued (and therefore invert¬ 
ible). Values of C D (() and c s (() in (I2.14dl) follow from a decomposition of the 
elastic strain energy 4 C(-)e:e into the deviatoric and the volumetric parts of the 
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strain tensor e. The stored energy of damage compliant with (12.14cl) is assumed in 
the form 

(4.2) b(C):=b iC, 

where b\ > 0 means the specific energy stored in the microcracks/microvoids created 
by damaging the material. By healing, this energy can be recovered back. The 
plastic yield stress compliant with (I2.14f|) is assumed in the form 

(4.3) er Y (C) = (ct Yi i-ct Y) 0 )C + cr Y ,o, 

where cr Yj i > a Y ,o > 0. The damage-dissipation potential is assumed in the piece- 
wise quadratic form 

(4-4) a(C) := ^(C+f + \a 2 {Cf + a 3 (C), 

where £ + = max{ 0 , £} and = max{—£, 0 } and 01 , 02,03 are given (material) 
nonnegative parameters. Values of ai and 02 determine rate-dependent parts of 
healing and damage model components and the value of 03 a rate-independent 
damage activation. The form of o(-) satisfies (12.14bl) . 

With respect to the fractional-step strategy of Section 3, we solve first for 
(u^. h , n^ h ) from the elastoplastic minimization problems (I3.5al) and then £* from 
the damage minimization problem (I3.5bl) recursively for k = 1, In view of 

the above shorcuts and simplifications, the minimization problems (13.5all and (13.5bl) 
rewrite as 

(4-5) fe, tt^) = argmin f f ^C^ 1 )(e(u+u* Th )-Tr) : (e(u+u* tTh )-n) 

- 9rh' U + CT Y (Cr/7 1 )l 7r - 7r rft 1 | N ) dx ~ f frh^dS, 

(4-6) Ch = argmin Q<C(C) (e(u^. h +u^ rh )-n^ h ) : (e(w^+w^ T/t -7r^)-6iC 

+ |«ivci 2 + ^a 1 (c-c\- 1 ) + + ^2(c-cvr+a 3 (c-c r Vr) 

where u is searched over Pl-elements satisfying Dirichlet boundary conditions, ir 
over PO-elements satisfying elementwise trace-free condition tr 7 r = 0 and £ over Pl- 
elements satisfying the nodal box constraint £ £ [0,1]. The form of (14.51) corresponds 
to the minimization problem of perfect plasticity with the elasticity tensor and the 
plastic yield stress depending on the damage variable in the previous time level. The 
energy in ( 1431 ) is transformed to an energy in the variable u only by substituting 
the elementwise dependency of 7 r on u , see BE] for more details. Then, the quasi- 
Newton iterative methods is applied to solve u^ h while n* h is reconstructed from it. 
More details on this specific elastoplasticity solver can be found e.g. in HEM- 
The damage minimization problem (14.61) represents a minimization of a nons¬ 
mooth but strictly convex functional. It can be reformulated to a modified problem 

(4.7a) argmin f ( ^-C(()(e(u* h +u^ Th )-Tr* h ) : (e(u^ h +u^ Th -n^ h ) 

C,z+,z_ Jn \ ^ 

- biC + —«|VC | 2 + ^a 1 {z + ) 2 + —a 2 (z _) 2 + a 3 Z- j dx, 
(4.7b) where z + = (C-^~ 1 ) + , Z-= 
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are additional ‘update’ variables. It should be noted that ( and ('^ 1 are Pl- 
functions and therefore Z+ and Z- are not P 1-functions in general on elements 
where nodal values of (—C,^ 1 alternate signs. However, if we restrict z + , Z- to Pl- 
functions while (14.7bl) is required on at nodal points, then (14.7al) actually represents 
a conventional quadratic-programming problem (QP), in which we require a linear 
and box constraints 

(4.8) C = Crh 1 + Z + ~ Z ~i Z + e [0> 1 - Crft 1 ]) z ~ e [0) Crh ^ 

A quadratic cost functional of this QP problem has a positive-semidefinite Jacobian, 
since there are no Dirichlet boundary conditions on the damage variable Q Note 
that the optimal pair (z+, Z-) must satisfy z+Z-=0 in all nodes, i.e. both variables 
cannot be positive. This can be easily seen by contradiction: If z+Z-> 0 in some 
node, then a different pair (z+— min{2L|_, Z-}, Z-— min{z + , Z-}) would again satisfy 
the constraints (14.811 but would provide a smaller energy value in (14.7al) . 

Our MATLAB implementation is available for download at Matlab Central as 
a package Continuum undergoing combined elasto-plasto-damage transformation , 
cf. [55]. It is based on an original elastoplasticity code related to multi-surface 
models [S]. The code is simplified to work with one surface variable only (which 
corresponds to the classical model of kinematic hardening) and sets the hardening 
parameter to zero to enforce perfect plasticity. It partially utilizes vectorization 
techniques of [43] and works reasonably fast also for finer rectangular meshes. 

5. Illustrative computational simulations. We consider a time-simulation 
of a 2-dimensional continuum visualized in Figure [5Tl describing two “plates” mov¬ 
ing horizontally in opposite directions with the constant velocity ±10 _8 m/s = 30 cm/yr. 
The model has applications in geophysics, specifically in modelling of tectonic and 
seismic processes in crustal parts of the earth lithosphere in the relatively short or 
very short time scales (meaning substantially less than a million of years) where 
the small-strain concept and solid mechanics are well relevant. The hardening is 
naturally considered zero. The damage variable is in the position of a so-called age¬ 
ing. The healing together with the damage-dependent plastic yield stress allow for 
periodically alternating fast damage and slow healing under external loading with 
constant velocity, which is a typical stick/slip-type events of flat partly damaged 
subdomains (so-called lithospheric faults) manifested by re-occurring earthquakes. 
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Fig. 5.1. Geometry used for the computational experiment, imitating the fault between two 
plates moving horizontally in opposite directions. The time-dependent Dirichlet conditions are 
prescribed on Td, using the constant velocity ±10“ 8 m/s = 30cm/year. 

The domain is assumed to be occupied by an elastic continuum specified 
by an isotropic homogeneous elasticity tensor in the form (14.11) with Ai = 7.5 GPa 
and fj,i = 11.25 GPa (which corresponds to Young’s modulus E Youn = 27 GPa and 
Poisson’ ratio v = 0.2 in the non-damage state) while the damaged material uses 
ten-times less moduli, i.e. Ao = 0.75 GPa and fio = 1.125 GPa in (14.11) . The yield 
stress cr y in (14.31) ranges between the values o Y , i = 2 MPa and o Y , o = ov,i x 10^ 12 . 
The damage-dissipation potential (14.41) is specified by constants ai = 100 GPa s 
and 03 = 10 Pa while the damage viscosity a 2 will vary. The stored energy of 
damage is b\ = 0.001 J/m 3 with the damage length-scale coefficient k = 0.001 J/m. 
The initial conditions ensure that ttq = 0 , Co = 1 (c* r Co = 1/2 in a middle narrow 
horizontal stripe). 
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The first numerical test is run for discrete times in the interval 0 < t < 400 ks 
with the equidistant time partition using the time-step r = 1 ks. The spatial dis¬ 
cretisation of the domain used a uniform triangular mesh with 4608 elements and 
2373 nodes; this mesh is available by setting ’level=2’ in the code (55], while finer 
uniform meshes can be generated by putting higher values of the ‘level’ parameter. 
Thus, 400 time-steps are computed and Figure [53] displays space-distributions of 
the shifted damage 1 — of the Frobenius norm of the plastic strain it, and of the 
von Mises stress | dev(cr)| at selected instants. 
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Fig. 5.2. Evolution of space-distributions of damage (the left column, displaying 1 —£), of the 
plastic strain (the middle column, displaying the Frobenius norm \ir\) and of the von Mieses stress 
(the right column, displaying |dev(ir)| / ). The displacement of the deformed domain is displayed 
magnified by the factor 12500. Distributions were computed for damage viscosity a 2 = 10 MPas. 

In order to see how the quality of discrete solutions depends on the time-step r, 
similar numerical tests are run for two additional time-steps r = 5 ks and r = 10 ks. 
The resulting energy balance (I3.9el) is displayed in Figure [531 Naturally, it is best 
fulfilled for the smallest considered time-step r = 1 ks. Figure 1531 shows the (hori¬ 
zontal component of the) reaction force which is here evaluated (very roughly) as an 
average from element values of von Mises stresses in the middle narrow horizontal 
stripe (i.e. the fault zone) shown in Figure [STTl A comparison of Figures [5731 and [531 
indicates that the energy balance (13.Dell is better satisfies in the purely elasto-plastic 
regime than within the undergoing damage. This becomes even more apparent if 
the damage process is speeded up by setting a smaller value ai = 0.1 MPas, cf. the 
left-hand parts of Figures 1531 and 1CT1 versus the right-hand parts. 

Dependence of the reaction-force evolution for varying viscosity of damage is 
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Fig. 5.3. Evolution of the stored and dissipated energy (=the left-hand side of (l3.9cl) for 
T varying) and the work of external loading (=the right-hand side of (13.9cD for T as a current 
time t) calculated for three different values of the time steps r = 10, 5, Iks, documenting the 
convergence of (l379cll towards the energy equality d2.Hct) proved in Proposition \3.3\ For less 
viscous damage this convergence is naturally slower than for a more viscous damage, cf. the left 
figure for a 2 = 0.1 MPas vs the right one for a 2 = 10MPas. 
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Fig. 5.4. Evolution of the reaction force corresponding to Figure 1 5. 3\ the time scales on the 
left and the right figures are different. Noteworthy, the force response is well converged even in 
situations when the energetics on Figure roi exhibits still big gaps. 


shown in Figure 15.51 for <22 as in Figures 15.3115.41 compared also with a smaller vis¬ 
cosity 02 = lkPas which already provides a response essentially identical to the 
even smaller viscosity 02 = 0.01 kPas (not displayed in Figure [575]) where conserva¬ 
tion of energy is numerically still more difficult to achieve. This indicates a certain 
tendency for convergence towards the model using rate-independent damage com¬ 
bined with rate-dependent healing (as in [371 Sect. 5.2.7]) and with perfect plasticity, 
which is theoretically not justified, however. 

Let us eventually remark that the a-posteriori information obtained from the 
residuum in the discrete energy balance (13.9cl) written at a current time t (as also 
used in Figure 15751) can be used to control adaptively the time step in a way to keep 
the numerical error in the energy under an a-priori prescribed tolerance and, on 
the other hand, not to waste computational time by making too small time steps 
in periods of slow evolution. We intentionally presented our numerical simulation 
on equidistant time partitions, but for actual geophysical simulations with very big 
difference in time scale between fast damage (earthquakes) and very slow healing, 
such an adaptivity is necessary. 
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Fig. 5.5. Dependence of the repulsive-force response on the viscosity of damage, the cases 
a ,2 = 10 and 0.1 MPas are (parts of) Figure \5J\ and are here compared also with even less viscous 
damage for 02 = 1 kPa s which gives essentially the same response as for the nearly inviscid 
case 02 = 0.01 kPas (not displayed, however); the time-step r = Iks. For decreasing viscosity, 
the rupture occurs earlier and propagates faster, showing a tendency to converge to an inviscid 
rate-independent (and theoretically not justified) damage model. 
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